Simulation of Droplet Trains in Microfluidic Networks 
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In this work we show that in a microfluidic network and in low Reynolds numbers a system can 
be irreversible because of hysteresis effects. The network, which is employed in our simulations, is 
taken from recent experiments. The network consists of one loop connected to input and output 
pipes. A train of droplets enter the system at a uniform rate, but they may leave it in different 
patterns, e.g. periodic or even chaotic. The out put pattern depends on the time interval among 
the incoming droplets as well as the network geometry and for some parameters the system is not 
reversible. 

PACS numbers: 



I. INTRODUCTION 

Fluids behave differently at the microscale. In this 
scale, factors such as surface tension, energy dissipation, 
and fluidic resistance are dominant and the flow stays 
laminar because of the small Reynolds number. These 
factors lead to unintuitive behavior in a microfluidic net- 
work which is a network of channels and pipes about 
one to hundreds of micrometer large jTJ [5] . Microfluidic 
devices have applications in inkjet printers, lab-on-a-chip 
devices, fuel-cells |3J, biochips |H[S] and recently, the pos- 
sibility of making logical gates and signal encoders and 
decoders has been studied [BJ [7J . 

In a recent work Fuerstman et al. look at the possi- 
bility of coding and decoding of signals using a simple 
network of microtubes which is made of just one single 
loop with two connecting pipes |7] (FIG.JT]). A constant 
flow of fluid enters the system from one side of the de- 
vice (input pipe) and leaves from the other end (output 
pipe) . This fluid is used as a medium to carry bubbles or 
droplets of an immiscible liquid. The droplets are large 
enough to block the pipes but small enough to be stable 
along the path. Reaching a T junction the bubbles must 
go only through one of the ways since they can not split. 
The decision depends on the pipes flux velocities which 
are history dependent and are related to the distribution 
of other droplets which have already arrived into the sys- 
tem. A droplet has the same cross section as the pipe, 
with a very small length, and in each cross section of the 
pipes, we have exclusively a droplet or the carrier fluid. 
Therefore, we can see droplets as binary signals which 
carry information on a suitable medium, and microflu- 
idic network can be represented by equvalent electrical 
circuits [5]. 

The fluid input flow in the system is in a constant rate, 
but the flow rate in any branch is a function of the loop 
geometry (length and diameter) and more interestingly, 
the presence of droplets in either branches. While the 
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capillary force of a droplet slows down the current in one 
branch, the other branch may become more favorable for 
the next coming droplet. Considering the input stream 
as a binary signal (droplet means 1, no droplets means 
0) droplets get different patterns (signals) in the output 
with respect to the time intervals between entering and 
the geometry of the loop. This can be seen as a way to 
code a signal. The reversible dynamics of the fluids at low 
Reynolds numbers suggests that the signal could be de- 
coded if the out put signal is fed to a system of pipes with 
the same geometry [TJ. The possibility of making logical 
microfluidic devices has also been examined by Prakash 
and Gershcnfcldof based on the same concept [6] . They 
have employed more complex networks to demonstrate 
how such a systems may act as logical gates. The possi- 
bility of microfluidic logic has raised also the interest of 
computer scientists [9]. 

Tracking the droplets inside the system in addition to 
time intervals between the output droplets, provides us 
information about the branch that a droplet selects. This 
gives us another informative binary signal which has been 
studied by Jousse et al. [TU] , 

We run a simulation based on the experiment [7J to in- 
vestigate the relation between the time intervals of the in- 
put signal and the pattern of the output. The simulation 
results however show that we are in low Reynolds num- 
ber regime which may show chaotic behavior because of 
the nonlinear history dependent dynamics of the system. 
Besides, the reversibility is also not a general property of 
the system and it may break even in non chaotic regime. 
We also see how related are the periodic patterns of time 
intervals and the path selection signals. 



II. THE MODEL 

Our microfluidic model network is simulating the ex- 
perimental network which has been introduced in refer- 
ence [7J. The flow enters the system through the input 
pipe. It reaches the first T junction and splits into two 
branches which join again at the second T junction. The 
two branches of the loop are identical except in their 
lengths. In fact, the length difference of the branches of 
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FIG. 1: The configuration of pipes in the experiment [7] and 
our simulations. From point A to C there are two ways, ABC 
and ADC. When the lenght of these two pathes are equal we 
call it a symmetric setup and when they are different, we 
name it asymmetric setup. 



the loop is one of our model free parameters (figure [I]). 

In each step of the simulation, we update the flow speed 
in all the pipes. As the geometry does not change dur- 
ing the simulation, the only relevant parameter which 
controls the flow velocity is the number of the droplets 
in the pipes. We can then simply follow the droplets in 
the pipes until the next upcoming event, in which one of 
the droplets reaches one of the bifurcation points. We 
recalculate the flux rates of the new configuration and 
continue in this way whether the droplet leaves or enter 
the loop. In a cylindrical pipe of length L and diame- 
ter D, at low Reynolds numbers, the flux Q is related to 
pressure difference between the pipe ends, AP, via 



AP 
~Q 



(i) 



where /i is the liquid viscosity. In the model, we also need 
to consider the effect of the droplets on the current. A 
droplet affects the flux velocity by introducing a capillary 
pressure drop, 



the branch that has a stronger flow and for the case of 
equal flux (in the range of our numerical accuracy), we 
suppose that one of the branches is always more favor- 
able. Then we recalculate the pressure differences and 
the flux rate after any event and continue in this way. 

In our simplified model, we study the effect of the 
length difference AL between the branches on the output 
signal. We also investigated the impact of time interval of 
input droplets, AT on the output signal, while all other 
system parameters are fixed. 

We perform the simulation for two different pipe set- 
tings simulating the original experiment. In one of them 
branches have the same length (AL — 0) with a ratio 
L/D = 40 (symmetric). In other one there is a differ- 
ence of 10% between the lengths with the same diame- 
ters (asymmetric). For both symmetric and asymmetric 
configuration we study the time coding of the output sig- 
nal as a function of the time intervals between entering 
droplets. 

To study the reversibility of the process we consider 
two different ways of signal decoding. In the both meth- 
ods the sequence of the time intervals between the output 
droplets is saved and used as input signal in the reversed 
process. In one of the methods which we call in-place 
decoding, we stop the coding process and reverse the di- 
rection of the movements. For instance, if the input flux 
is entering from the left and exiting from the right, we 
reverse the flux direction from right to left and the distri- 
bution of the droplets inside the system also is preserved 
at the time of the flux reversion. In the other decoding 
method we suppose that we have only the coded signal 
and a decoding device, identical to coding device. That 
means we feed the coded signal to an initially empty sys- 
tem with exactly the same geometry of the coding device 
but in the reverse order. From the application viewpoint, 
this method may be more interesting and in contrast with 
the first method we call it out-of-place decoding. 
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are surface tension between the gas and the liquid and 
capillary number, respectively [TTJ Q2] ■ Here v indi- 
cates the droplet velocity which is the same as that of 
the carrier fluid. 

According to the above equations, we are able to calcu- 
late the pressure difference in any branch. Flux conser- 
vation on the network nodes results in a system of cubic 
equations for the flux rates which could be solved numer- 
ically. The flux rates of the pipes are preserved as long 
as the number of droplets in all the pipes are unchanged. 
Therefore, we calculate the position of the droplets over 
the time until a new droplet enters a pipe or a moving 
droplets meet the end of a pipe (we call it a new event) . 
When a droplet reaches the bifurcation point, it selects 



III. THE RESULTS 

In the simulation we enter the droplets in an uniform 
flat rate. Then an input signal is simply specified by a 
single parameter AX^, the time interval between entering 
droplets. But in the output the time intervals between 
the droplets might be different. Then the output signal 
is specified by a sequence of AT n = T n — T„_i where n 
counts the droplets. 

For large enough ATi's it is trivial that the output pat- 
tern is identical to the input pattern. By decreasing ATi, 
at a certain point, we start to get different patterns at the 
output. This happens if a droplet reaches the entering T 
junction, while the former droplet has not left the loop 
yet. So the droplet at the junction will select the other 
branch which is empty and has a higher flow rate. At this 
point we start to get a periodic output pattern (periodic 
sequence of AT n s) that are different from the uniform 
input signal. The period length (the number of droplets 
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FIG. 2: Output time intervals versus input time interval 
(top) for symmetric setup (left), with two identical branches, 
and asymmetric setup (right), in which, one branch is longer. 
The period of the output signal as a function of the input 
time interval is shown below each case. 



in a periods) increases as we decrease ATj. By even more 
decreasing ATi for some smaller time intervals, we start 
to get disordered patterns with no visible order or peri- 
odicity on the output (FIG{2]), there, the period length 
reaches the number of droplets involved in the system 
which means the system is either aperiodic (chaotic) or 
has a period longer than the simulation time. 

Similar to many other chaotic systems, here we also ob- 
serve that the system reemerge into the chaotic regime in 
different points for the both symmetric and asymmetric 
models. Interestingly, there are some intermediate inter- 
vals where, the coder loop is invisible and the output and 
input are identical. The period of the output signal as 
a function of input time interval is shown in FIG |2] and 
FIGJH 

As it is shown in figure Q, similar behavior has been 
observed in both symmetric and asymmetric networks 
however the figures are different. This is indicating that 
even at low Reynolds numbers, the system shows chaotic 
behavior merely as a result of history dependent dynam- 
ics of the model. The same feature have been observed 
in history dependent ecological models |14j . 

Similar to many other nonlinear systems, the output 
signal may change in time and we need to wait until it 
converge to its attractor. To sample the output intervals, 
we wait until a stationary distribution function of the 
intervals is achieved. This might take long time. 

We study the reversibility (decoding) of the process in 
both in-place and out-of-place methods. It is not surpris- 
ing that in the chaotic regimes, the process is decodable 
with neither methods. It shows that in the chaotic regime 
the input signal information is lost and it can not be de- 
coded. Instead it is possible to restore the original signal 



when we have periodic outputs however it is not always 
true that all the processes are reversible when they are 
not chaotic. This is more surprising for the case of in- 
place decoding. 

The different behaviors of the system for some input 
time intervals are presented in figure [3] for asymmetric 
setup. In the figure the patterns of time intervals are 
shown for "before" and "after" the reversion of the flow. 
First, the flow goes normally through the network and 
the output is shown for some droplets (the left part of 
each graph), then we look at the time interval between 
these droplets when the flow is reversed (the right side 
of each graph). The graphs differ only in the time in- 
tervals of the input signals. Therefore, as it is shown in 
the figure, even for periodic output signals, the initial 
signal sometimes is reconstructable and some times not. 
It is also possible that a periodic encoded signal shows 
chaotic behavior in its return (figure |3|;). The similar 
behavior is observed for symmetric setup also. Here as 
we are at low Reynolds numbers, time reversality might 
be expected. But again history dependence of the pro- 
cess plays its role. A droplet chooses a branch that has 
a faster flow. But when it leaves the branch, the flows 
might have been changed. As a result, when we re vers 
the time, the droplet who has just left a branch, might 
decide to go in the other branch for it has a faster flow. 

Looking at figure|3]one might ask if the decoded signals 
are stationary solution or not. This figures have been se- 
lected for demonstration because of fast convergence of 
solution, however, this is a relevant question. Thus we 
do in-place decoding for the whole range of ATi with 
appropriate witting time both in the encoding and the 
decoding process. The map of encoded signals for both 
symmetric and asymmetric models have been shown al- 
ready in figure [2] The corresponding decoded signals are 
shown in figure 0] For both models there are some area 
that the signals are decodeable. 

In the out-of-place decoding method, the signal can 
again be decoded for some periodic signals. As we have 
a periodic signal, there is an additional degree of freedom, 
the phase. That means when we are going to feed a pe- 
riodic signal into an empty device, which droplet should 
go first. We observe that for the periodic signals the re- 
sult is sensitive to this phase (the index of the starting 
droplet). Then in a periodic pattern, if we start decod- 
ing from different phases, we might get different signals. 
The resulting signals are again periodic and the original 
signal may be among them. 

An example is presented in figure [5] The figure [5^ 
shows a periodic signal which is encoded by asymmetric 
setup with a flow of droplets with ATi — 1/15. The num- 
bers in the figure indicate the phase. As the signal has a 
period of 5, we can feed this sequence into an empty de- 
coding device starting from any of this droplets. As the 
figure shows different patterns are archived if we start 
decoding process with droplets numbered 1 and 4, while 
the original sequence is constructed if we start from the 
other three droplets. The similar futures are observed for 
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FIG. 3: In-place decoding for asymmetric model: the time 
intervals are shown for "before" and "after" the reversion of 
the flow for 300 (a and b), 20000 (c) and 3000 (d) droplets, 
in the left and the right side of each graph (see text). For 
AT — 1/10.00 (a), there is a regular and periodic coded 
pattern which is decoded perfectly. For AT = 1/11.45 (b), 
there is again periodic coded output, but the time reversal 
pattern is a different periodic signal. For ATi = 1/42.60 (c), 
however, the coded signal is periodic, in reversion there is 
no attractive pattern and the system shows chaotic behavior 
very soon. For AT — 1/54.00 (d and e), the output pattern 
is aperiodic with the values concentrated on some strips. The 
reversion produces a similar pattern with more fluctuations. 



symmetric model, however, only the results of asymmet- 
ric are presented here. 

Up to now we only studied the output signal. We can 
go furthere and look at the intermediate pathway and 
study the patterns of pipe selections. We can follow the 
train of the droplets in input we can follow them and see 
any droplet is choosing which branch, say up (t) or down 
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FIG. 4: Time intervals between the droplets of in-placed de- 
coded signals as a function of uniform time intervals between 
the droplets in original signal (top) and the period of the out- 
put signal as a function of the input time interval (bottom). . 
For both symmetric (left) and asymmetric (right) setups there 
are some regions where the signal is decoded successfully. The 
corresponding encoded signals are shown in figure [2] 



(|). This give us another way of decoding the signal. 
This definition of signal has been introduced by Jousse 
et al. [10 in their experiment. We shall try the same with 
our model and compare the results with the experiment. 

Here we look at the signal in our asymmetric model 
and without lack of generality, we suppose the shorter 
branch is labeled |. If we count the number of droplets 
passing through "up" and "down" branches and show 
them by and n^, respectively, then p = n.|/(n.| + njj 
and q — ri|/(ri| + njj are probability of a bobble pass- 
ing through "up" or "down" branches. We expect in a 
long run, more droplets pass through the shorter branch 
(|) than longer (j.), then I > r = p — q > 0. But r is 
not only related to the pipe lengths, but also it varies by 
ATi. For a large enough ATi, when any droplet reaches 
an empty loop, all droplets pass through shorter branch 
and r = 1. Decreasing ATi in some point the droplets 
choose "up" and "down" branches alternatively. Here we 
have r = 0. One should note that, here the up-down sig- 
nal is alternative and sequence of output intervals is still 
uniform and identical to input signal. Thus if up-down 
signal is periodic the sequence of time intervals is also 
periodic, but they don't need to have the same period. 
As we have considered very small difference between the 
pipe lengths, for any periodic up-down signal and in any 
period the difference between the numbers of j and J, is 
or 1. That means for any periodic signal with period 
P, r is either or 1/P, for even and odd P numbers, 
respectively. 
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IV. CONCLUSION 



In this work we used a simple setup of microfiuidic 
network (FIG. [IJ and a set of deterministic equations to 
simulate the passage of droplets or bubbles in the mi- 
crofiuidic network. The fluid had a given constant flow. 
Which carries droplets. The time interval between suc- 
cessive droplets is constant at the input. We found that 
patterns of droplets at the output depends on the time 
interval between the input droplets. We observed that 
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FIG. 6: Shows the ratio of the difference between the number 
of droplets passed through each pipe versus input time inter- 
val (no — ni)/ ( n o + Tii). Symmetric setup (left), Asymmetric 
setup (right). 







10 15 20 25 



2.5 
2 

^1.5 

C 

% 1 

0.5 

0, 




10 15 20 25 



FIG. 5: Out-of-place decoding: a periodic output signal (top) 
of asymmetric setup for ATi = 1/15 has a period of 5. The 
output depends on the phase of the signal in decoding process, 
if we start from droplets 1 (middle) or from for 4 (bottom). If 
we start with the others (the numbers 2, 3 or 5) the original 
uniform signal is reconstructed 



in general we have both periodic and chaotic patterns of 
output depending on the rate of entering droplets. 

We found that although deterministic the system is not 
always reversible. We reversed the direction of time (or 
flow) in simulations and we found that in some situation 
the process is irreversible and droplets choose the other 
branch to come back which means a droplet leaving one 
of the branches, may choose the other one to return, if 
the time is reversed. This also depends on the input time 
interval of droplets. 

We observed that there is a relation between irre- 
versibility and chaos in our system. In general, when 
the output signal is periodic, it may be decoded to the 
original input signal by just reversing the flow direction. 
When it is chaotic, it is not possible to generate the input 
signal (decode) . 
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